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Abstract 

Nuclear structure quantum Monte Carlo methods such as Green's function or auxiliary field 
diffusion Monte Carlo have used phenomenological local real-space potentials containing as few 
derivatives as possible, such as the Argonne-Urbana family of interactions, to make sampling 
simple and efficient. Basis set methods such as no-core shell model and coupled-cluster techniques 
typically use softer non-local potentials because of their more rapid convergence with basis set size. 
These non-local potentials are usually defined in momentum space and are often based on effective 
field theory. Comparisons of the results of the two types of methods can be difficult when different 
potentials are used. We show methods for evaluating the real-space imaginary-time propagators 
needed to perform quantum Monte Carlo calculations using such non-local potentials. We explore 
the universality of the large imaginary time propagators for different potentials and discuss how 
non-local potentials can be used in quantum Monte Carlo calculations. 
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I. INTRODUCTION 



Quantum chromodynamics (QCD) is the fundamental theory that underlies the descrip- 
tion of atomic nuclei, where quarks and gluons are the degrees of freedom. However, despite 
four decades of studying QCD, little connection has been made between QCD and low-energy 
many-body nuclear dynamics. The most direct computational application of QCD, lattice 
QCD, has made much progress in the past decades but remains some distance away from 
being a practical tool for computing many-body nuclear observables. For recent reviews of 
the outlook for lattice QCD calculations as they apply to nuclear physics see Refs. 
While direct applications of effective field theory have made progress [3H6| the characteriza- 
tion of atomic nuclei in terms of phenomenological two- and three-body nucleon interactions 
remains the standard starting point for most nuclear structure calculations. 

Quantum Monte Carlo is one of the most successful methods for nuclear matter and 
nuclear structure calculations. Green's function Monte Carlo has solved for many low lying 

n n 

states of nuclei for A < 12 [7J, |8[ and auxiliary-field diffusion Monte Carlo can calculate 
much larger nuclei and nuclear and neutron matter js-ll]. These methods have used phe- 
nomenological local real-space potentials containing as few derivatives as possible, such as 



the Argonne-Urbana family of interactions 



Ref. \S 



12Ml4i] . to make sampling easy and efficient. See 



for a review of Green's function Monte Carlo results. 



However, there are other successful approaches that can reach A = 12 and beyond. Basis 



set methods such as the no-core shell model 



15 



, ll6| and coupled-cluster techniques 



17- 



have typically used softer non-local potentials as these have more rapid convergence with 
basis-set size. These potentials are typically defined in momentum space and are often 
derived from chiral effective field theory such as the next-to-next-to-next-to-leading order 
(N 3 LO) interaction of Ref. Q. 

Comparison of the results of basis set methods with Monte Carlo calculations can be 
difficult when different potentials are used. Quantum Monte Carlo methods use propaga- 
tion in imaginary time to project out the low-energy states of a quantum many-particle 
system. The propagation is performed by first writing the many-body short-imaginary-time 
propagator. Accurate methods use a pair-product approximation 2l|, |22j where the many- 
body propagator is written in terms of the propagator for each pair of particles. In order to 
perform a quantum Monte Carlo calculation using an arbitrary pair potential, we need to 



calculate the pair propagator in imaginary time. 

The aim of this paper is to show how to evaluate this real-space imaginary-time pair 
propagator needed to perform quantum Monte Carlo calculations using non-local potentials 
(Sec. II), to demonstrate the consistency of the propagators using different potentials at 
large imaginary times (Sec. Ill), and to discuss how to use these propagators to calculate 
the properties of nuclei and nuclear matter using these non-local potentials with quantum 
Monte Carlo methods (Sec. IV). 



II. METHODS 



Non-local potentials are typically generated in momentum space from an effective field 
theory. With a potential defined in momentum space, V(k,k'), it is natural that we pro- 
ceed by constructing a Hamiltonian in momentum space as well. Our normalization and 
completeness conventions for our continuous real-space and momentum-space basis states 
are 

1 = J dh\r){r\ = J ^\k)(k\. (1) 

r is the separation vector of the two nucleons and p = hk the conjugate momentum. The 
overlaps between the states are 

(r\k) = e lkr . (2) 

We work in the standard channel basis where J 2 , J z , L 2 , S 2 , S z , T 2 and T z are good quantum 
numbers, with the total spin, S = Si + S 2 , the total angular momentum J = L + S, and 
total isospin T = T x + T 2 . We choose our basis states as \r.JMLSTT z ) and \k.JMLSTT z ) , 
with normalization and completeness given by (suppressing J, S, T, and T z ) 

1 = Yl / r 2 dr\rLM)(rLM\ 

=y, r ^ kLM){kLMi (3) 

LM J ° ^ > 

The overlaps are 

(rLM\kL'M'} = 4m L j L (kr)5 LL ,5 MM , (4) 

For numerical work, we compactify our real space to a sphere of radius R. We choose 
the Dirichlet boundary condition on the sphere which forces our momentum-space spectrum 
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to be discrete with k^R being the zeros of the spherical Bessel functions, j^k^ R) = 0. 
Below we often drop the superscript (L) when its value is clear from context. The discrete 
momentum states \k n LM) are chosen with unit normalization so that 



(rLM\h n L'M f ) = ^ — ^ R ^ 3dk n r)5 LL ,5 M w ■ (5) 

Our transformations can now be treated as orthogonal-matrix multiplications. 

The momentum-space Hamiltonian for the uncoupled channels where L = J (for a given 
set of the quantum numbers — J, M, L, S, T, and T z — which we suppress below) is 

h 2 k 2 

(k m \H\k n ) = — 5 mn + V(k m , k n ), (6) 
lm r 

with m r the reduced mass. For the coupled channels where the potentials couple the L = 
J ± 1 states together the Hamiltonian is 

(k m L\H\k n L>) = ( 2 - 5 ™ ' ^ } ^ + /- +( km ' (+) , (7) 

\ V-\ — (km , k n ) 2 m r — ^mn + ^++(^m , k n ) I 

where the superscripts (— ) and (+) correspond to L or V having values of J — 1 and J + 1. 
We then construct the momentum-space, imaginary-time propagator by diagonalization of 
the Hamiltonian, giving 

(k m \e~ Hr \k n ) = ^(U^>e~ £lT (^IU, (8) 
i=i 

for the uncoupled channels, and 

(k m L\e- Hr \k n L') = J2(kmL\i>i)e- EiT (il>i\knL'), (9) 
i=i 

for the coupled channels. The {\ipi) } are eigenvectors of the Hamiltonian with corresponding 
eigenvalues {Ei}. Nk is the number of discrete momentum states we keep. We ensure that 
N k is large enough such that the propagators converge. An estimate of how large k max 
should be can be given by considering the kinetic energy alone. We want A; max such that 
exp (^— h 2^ a * r ) can be neglected. In practice, we check the convergence by doubling our 
estimate for k nmx and ensuring our results do not change to the desired precision. With 
these methods, it is easy to ensure that the numerical truncation errors are completely 



negligible. For example, for the results shown here, we use k max = 40 far 1 (N k ~ 80), and 
the propagators have truncation errors less than 10~ 10 . 

After transforming to real space, we have the matrix elements 

(r JMLSTT Z \e~ Hr \ r' J M L' STT Z ) . (10) 

However, for use in Monte Carlo codes, we want the propagators in a 3D real-space basis, 
\r9<pSM s TT z ), 

\r9(j>SM s TT z ) = cJ sm s lmJlm l [QA)V JMLSTTZ . (11) 

JMLM L 

C is a Clebsch-Gordan coefficient, Y a spherical harmonic. 

In quantum Monte Carlo calculations, the particle positions are typically sampled from 
the central part of the propagator. The non-central parts are then included in the spin- 
isospin samples (auxiliary field diffusion Monte Carlo) or sums (Green's function Monte 
Carlo). We will sample the propagators for these non-local potentials in the same way. 
Here, we define the central part of the propagator as the trace over all spins and isospins. 
For convenience we also choose a particular coordinate system where the initial separation 
lies along the z axis, and the final separation is in the xz plane such that we may take 
6 = <j) = 0' = and we can visualize the central part of the propagator as a function of 
r — r' and 9'. For any particular application, we can always rotate into this configuration, 
propagate, and rotate back. The central part of the propagator is then written as 

G(r,r',0';r) = {rSM s TT z \e' HT \r'd' SM S TT Z ). (12) 

SM S TT Z 

III. CONSISTENCY OF THE PROPAGATORS AT LARGE IMAGINARY TIMES 

In the limit of large imaginary times, we expect that the propagators for different po- 
tentials should agree. The propagators are essentially density matrices for the two nucleon 
system: 



P = „-e-t ' tr P = ( 13 ) 



corresponding to thermal equilibrium at the temperature k B T = r _1 . Now, since any 
measurable quantity can be written as an expectation of a Hermitian operator O which 
can be obtained via (O) = tr(pO), the density matrices (propagators) we obtain for the 



various potentials contain all the measurable information for this system. If the position or 
momentum of the nucleons could be determined with arbitrary precision, the density matrix 
would be in principle measurable. Since the position and momentum are not well defined 
for arbitrary values, the propagator is not completely measurable. If the various potentials 
we use are phase-shift equivalent — meaning they reproduce the physical scattering data 
at or below £?i a b ~ 350 MeV (£ c . m . ~ 175 MeV) — then we would expect that starting at 
imaginary times r « (175 MeV) -1 , the various propagators should begin to agree more and 
more. In fact, we find that for r rj (50 MeV) -1 , the higher-energy modes not constrained 
by current experimental data do not contribute substantially to the propagator. 

We see from Eq. ( I13p that at r ^ (175 MeV) -1 , energies of 175 MeV and above are 
suppressed by a factor of 1/e. Therefore, it is not surprising that at r w (50 MeV) -1 , 
where energies of 175 MeV and above are suppressed by a factor of 1/e 3 ~ 0.0498 we find 
relatively good agreement between the propagators with different potentials. This result is 
analogous to the renormalization group results leading to the V\ ow & potential of Ref. [23] 
and the similarity renormalization group results of Ref. j^j, where the high energy modes 
are integrated out. 



Figures [THH] demonstrate these findings for three potentials: Argonne Vx& (AV 18 ) [12], 
N 3 LO, and N 3 LO(600) [^(J. For the diagonal cases, where we take k = k', we define a 
quantum potential, V q (k,k;r) through the equation 

(k\e~ HT \k) = (k\e- Ho $e- v « {k ' k - T)T e- H ^\k), (14) 

with Hq the free-particle Hamiltonian 

"o = T=^. (15) 

In the off-diagonal cases (k ^ k') we choose a particular k value and plot against k! . For 
visual comparison, we subtract the free-particle propagator, g(k, k') — go(k, k'), since at the 
point k = k', the kinetic energy component is large and obscures the result. 

Figures HHH show the quantum potential in the singlet ( 1 S'o), uncoupled triplet ( 3 -Po); 
and coupled triplet ( 3 Si and 3 Di) channels. The imaginary times chosen correspond to a 
typical time step used in Green's function Monte Carlo calculations, r = (2000 MeV) -1 and 
imaginary times that roughly correspond to center-of-mass energies of 350 MeV, 175 MeV, 
and 50 MeV. As we have discussed above, the imaginary time of r = (50 MeV) -1 is the time 
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at which the effects attributable to energies of 175 MeV and above are effectively integrated 
out. It is interesting to note that the agreement of the quantum potentials is only good 
up to k ~ 2 fm" 1 : this is the approximate momentum value, k, one would associate with 
the corresponding kinetic energy: 4^ = 175 MeV. This relationship (better and better 
agreement — but only up to some cut off — as the potential is evolved) is precisely what is 
found in similarity renormalization group and V\ ow & approaches. 

It is tempting to interpret r as an evolution parameter for the quantum potential in 
the same sense that the similarity renormalization group approach has s or A (see, for 



example, 24j])- However, it is not clear if a direct comparison is at all trivial. In the 
similarity renormalization group approach, the evolution parameter s (and therefore A, since 
A = l/s 1//4 ) is defined through the rather simple evolution equation for the potential 

where G s is a Hermitian operator that generates the transformation. (G s is often chosen 
to be T, the kinetic energy). If we view our quantum potential akin to the similarity 
renormalization group's V s , our defining equation is 

e -Hr = e - T r e -V q (r)r e -TZ 

We can expand this using the Baker-Campbell-Hausdorff relation which gives an infinite 
series of nested commutators. The lowest order terms in an expansion in r are 

V q (r) =V-& ([V, [V, T)) - i[T, [T, V))) + ■■■ (18) 



where the double commutators are suggestive of Eq. (Tl6l) . but clearly not the same. The 
differential equation satisfied by V q {r) is not a simple, compact expression. 

Figures |5]|9l show the off-diagonal elements of the singlet ( 1 So) ) uncoupled triplet, ( 3 Po), 
and coupled triplet ( 3 Si, 3 D ly and 3 Si- 3 Di) channel propagators. As discussed above, they are 
shifted by the free-particle result to make comparisons easier. What we can see from these 
figures is a general trend towards universality with at least two caveats. First, the propa- 
gators tend to converge most rapidly around the point where k = k'. Our interpretation of 
this result is that low momentum transfer behavior is constrained by the phase-shift equiv- 
alence of the various potentials whereas higher momentum transfer behavior is not. Second, 
some channels converge better than others. For example, Fig. [8] has still not converged 



at r = (50 MeV) -1 , and indeed, does not appear to converge well until r = (10 MeV) -1 . 
This may be due to this channel's sensitivity to the tensor part of the interaction which 
the different potentials treat differently. These differences may point to true, quantifiable 
distinctions between the potentials. 

Even though the potentials give propagators with the same sort of long imaginary time 
behavior, the many-body physics of the nucleus may not allow the use of the propagators 
in this regime. As mentioned above, in Green's function Monte Carlo calculations, the 
imaginary time step needed to accurately approximate the many-body Green's function by 
the pair-product is r = (2000 MeV) -1 . For larger time steps, commutator terms in the 
Trotter breakup spoil the approximation. Since, in the pair product approximation, these 
commutators occur only when three nucleons are close together, this indicates that, three- 
body effects will also be important. That is, to be able to use the propagators in the limit 
where they become model independent, would likely require that three- and more-body 
terms in both the interaction and the propagators be included. This likely means that 
use of potentials like V% w & for many-nucleon calculations will need to include many-body 
interactions. 



IV. QUANTUM MONTE CARLO WITH NON-LOCAL POTENTIALS 

We now turn to the central parts of the propagators for the non-local N 3 LO and 
N 3 LO(600) potentials we have been considering throughout. Since the central part of the 
propagator is sampled in a quantum Monte Carlo calculation it should be positive-definite 
to avoid sign problems. 

Figures [TU1 and [TT1 plot the central part of the propagator in real space, where the coordi- 
nates x' and z' are such that the origin corresponds to the final separation equal to the initial 
separation. That is, the relative coordinates are equal: r = r'. The precise transformation 
between the original coordinates, r, r', and 9' and the new coordinates, x' and z' is given by 

r ' 2 = x ' 2 + ( r + z ') 2 (19a) 

COS 0< = ^ = J + Z ' (19b) 
^x' 2 + (r + z') 2 

and can be visualized in Fig. [12j 

The structure we can see is a Gaussian-like peak about the initial separation as well as 
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an antisymmetric trough at the position that corresponds to the two nucleons undergoing 
a position interchange: r' = — r. The antisymmetric point is built in from tracing over the 
spins and isospins as in Eq. (I12p . and would be present even for Argonne t>ig. This point 
gives no extra difficulty — since it comes from the fermion character of the nucleons, it 
will be dealt with in the same way that the fermion sign problem is dealt with in Green's 



function or auxiliary field diffusion Monte Carlo. That is, a path constraint 



25j is imposed 



that eliminates the fermion sign problem. The constraint can then be released and forward 



walking steps taken [25l . |26| to improve the results and check the effect of the constraint. 

However, if we zoom in on the shifted origin, as in Figs. [121 and dU setting any positive 
parts of the propagator to zero, and make sure we are clear of the antisymmetric region, 
we find that the propagators appear to be "ringing", much like Friedel oscillations. These 
negative parts may make it more difficult to perform quantum Monte Carlo calculations and 
keep the sign problem under control. However, these negative parts are quite small, of order 
1CT 1 fm -3 , whereas the peak of the propagator is of order 10 2 fm -3 . In fact, a typical slice 
through the propagator in the x' direction looks like Figs. [15] and [16j In many cases, the 
negative parts are negligible. 

A straightforward method to take the small negative regions into account, is to first set 
any negative part (not associated with the antisymmetry) to zero, run a quantum Monte 
Carlo calculation until it converges, and then add the negative parts back in, in a perturba- 
tive fashion, using forward walking exactly as for the fermion sign problem. The extra sign 
changes from the propagator are handled in the same way as sign changes from the fermion 
character. 

We can estimate the fraction of walkers in the initial time step that may be given negative 
weights by comparing the integral of the absolute value of the propagators to the integral 
of the propagators. That is, we estimate the fraction of walkers with negative weights, / by 

, = j<Pr'\[\G{ry-r)\-G{ry-r)} 

1 Jd 3 r'\G(r,r';r)\ ' 1 ' 

We calculate the integral over the upper half volume to exclude the interchange. For N 3 LO 
with an initial separation of 1.0 fm, we find / ~ 0(1O~ 2 ). If we now take the very conserva- 
tive estimate for the alpha particle that all six pairs may be this close at one time and that 
the negative weights are acceptable so long as the fraction of walkers with negative weights 
is less than 1/e, we find we can take approximately ten steps for N 3 LO [with time step 



(2000 MeV)" 1 ]. Green's function Monte Carlo typically uses forward walking of about 10 
to 20 steps of (2000 MeV) -1 25|, |26|. Therefore forward walking will allow us to remove any 
bias from the negative parts of the propagator. This can be compared to forward walking 
keeping the propagator constraint, but releasing the fermion constraint to separate the two 
effects. 

In the above analysis, we have assumed that the imaginary time step used will be the 
same as that used in Green's function Monte Carlo calculations with the Argonne family 
of potentials. Since the N 3 LO potentials are softer, the relevant commutators terms will be 
smaller, and longer time steps may be possible. For longer time steps there is much less 
ringing, and the calculations will be substantially easier. 



V. CONCLUSIONS 

We have shown how to calculate the imaginary time pair propagators needed for quantum 
Monte Carlo calculations of nuclei and nuclear matter using non-local potentials in momen- 
tum space. The method is general enough to handle any non-local potential in momentum 
or real space, but in this paper, we focus on those derived from effective field theory (N 3 LO). 

We find that the propagators display a universal behavior at large imaginary times, 
consistent with our expectations from renormalization group methods and the fact that the 
potentials are phase-shift equivalent, meaning that they reproduce the scattering data at or 
below laboratory energies of 350 MeV. 

The central propagators sampled during Monte Carlo simulations for local potentials 
are expected to be positive-definite. Without this property, sign problems can develop. 
We find that for N 3 LO with a 500 MeV or 600 MeV cutoff in momentum space, the central 
propagator is not positive definite. However, the negative parts consist of "rings" reminiscent 
of Friedel oscillations and their magnitude is quite small compared with the overall shape 
of the central propagator. Since these potentials were developed in momentum space, no 
attempt was made to influence their behavior in position space. It may be possible to modify 
the N 3 LO potentials in such a way that they continue to reproduce the Nijmegen data with 
a low x 2 ; are still relatively soft, but have reduced ringing behavior. A modification of the 
choice of the regulator function used in the calculation of the N 3 LO potentials: V(k, k') — > 
V(k, k')e~( k / A " >2 '' e~( k ' ' l A -) 2v ; where A is the cutoff value, and v is the order of the calculation, 
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[y = 4 for N 3 LO) may help. In any case, quantum Monte Carlo calculations should still be 
possible by using a modified path constraint as described above, and we are implementing 
these calculations. 

While we have concentrated on calculating the imaginary time pair propagators from 
phenomenological potentials, it is amusing to note that since the few-body imaginary time 
propagators are simply imaginary-time correlations of the appropriate nucleon operators, 
they might eventually be directly extracted from lattice QCD calculations. 
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FIG. 1. (color online) The quantum potential in the So partial wave in momentum space as the 
propagator is calculated for successively longer imaginary times. 



13 



4 



x=(2000 MeV)" 1 


Al/ 18 — 
N 3 LO 
N 3 LO(600) 




x=(350 MeV)" 1 


i ■ i ■ i 

x=(175 MeV) " 1 


, 1 , 1 1 1 1 

x=(50 MeV)" 1 



1 2 3 1 2 3 4 

/c(fm" 1 ) /((fm" 1 ) 

FIG. 2. (color online) The quantum potential in the 3 Po partial wave in momentum space as the 
propagator is calculated for successively longer imaginary times. 
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FIG. 3. (color online) The quantum potential in the 3 S\ partial wave as the propagator is calculated 
for successively longer imaginary times. 
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FIG. 4. (color online) The quantum potential in the 3 D\ partial wave as the propagator is calculated 
for successively longer imaginary times. 
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FIG. 5. (color online) The off-diagonal momentum-space propagator (minus the free-particle prop- 
agator) in the partial wave as the propagator is calculated for successively longer imaginary 
times. 
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FIG. 6. (color online) The off-diagonal momentum-space propagator (minus the free-particle prop- 
agator) in the 3 Po partial wave as the propagator is calculated for successively longer imaginary 
times. 
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FIG. 7. (color online) The off-diagonal momentum-space propagator (minus the free-particle prop- 
agator) in the 3 Si partial wave as the propagator is calculated for successively longer imaginary 
times. 
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FIG. 8. (color online) The off-diagonal momentum-space propagator (minus the free-particle prop- 
agator) in the 3 Di partial wave as the propagator is calculated for successively longer imaginary 
times. 
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FIG. 9. (color online) The off-diagonal momentum-space propagator (minus the free-particle prop- 
agator) in the s Si- 3 Di partial wave as the propagator is calculated for successively longer imaginary 
times. 
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FIG. 10. (color online) The N 3 LO central propagator for initial separation, r = 0.5 fm, and 
imaginary time, r = (2000 MeV) -1 . 
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FIG. 11. (color online) The N 3 LO(600) central propagator for initial separation, r = 0.5 fm, and 
imaginary time, r = (2000 MeV) -1 . 
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FIG. 12. The coordinate system used to visualize the central propagators. O is the original origin, 
and O' the shifted origin such that r = r' corresponds to the shifted origin. 




FIG. 13. (color online) The negative parts of the N 3 LO central propagator for initial separation, 
r = 0.5 fm, and imaginary time, r = (2000 MeV) -1 . 
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FIG. 14. (color online) The negative parts of the N 3 LO(600) central propagator for initial separa- 
tion, r = 0.5 fm, and imaginary time, r = (2000 MeV) -1 . 
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FIG. 15. A typical slice through the N 3 LO central propagator for initial separation, r = 0.5 fm, 
and imaginary time, r = (2000 MeV) -1 . 
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FIG. 16. A typical slice through the N 3 LO(600) central propagator for initial separation, r = 0.5 
fm, and imaginary time, r = (2000 MeV) -1 . 



25 



